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We investigate the existence and properties of equipotential surfaces and La- 
grangian points in non-synchronous, eccentric binary star and planetary systems 
under the assumption of quasi-static equilibrium. We adopt a binary potential 
that accounts for non-synchronous rotation and eccentric orbits, and calculate 
the positions of the Lagrangian points as functions of the mass ratio, the degree 
of asynchronism, the orbital eccentricity, and the position of the stars or planets 
in their relative orbit. We find that the geometry of the equipotential surfaces 
may facilitate non-conservative mass transfer in non-synchronous, eccentric bi- 
nary star and planetary systems, especially if the component stars or planets are 
rotating super-synchronously at the periastron of their relative orbit. We also 
^ ■ calculate the volume-equivalent radius of the Roche lobe as a function of the 

four parameters mentioned above. Contrary to common pract ice, we find that 



replacing the radius of a circular orbit in the fitting formula of lEggletonl (119831 ) 
with the instantaneous distance between the components of eccentric binary or 
planetary systems does not always lead to a good approximation to the volume- 
equivalent radius of the Roche-lobe. We therefore provide generalized analytic 
fitting formulae for the volume-equivalent Roche lobe radius appropriate for non- 
synchronous, eccentric binary star and planetary systems. These formulae are 
accurate to better than 1% throughout the relevant 2-dimensional parameter 
space that covers a dynamic range of 16 and 6 orders of magnitude in the two 
dimensions. 



Subject headings: Celestial Mechanics, Stars: Binaries: Close, Stars: Planetary 
Systems 
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Introduction 



The Roche model has served for a long time as a fundamental tool in the study of the 
interactions and observational characteristics of the components of gravitational two-body 
systems. Among the most noteworthy applications are the modeling of the shapes of binary 
components, the study of mass transfer in interacting binaries, and the orbital stability 
of satellites in planetary systems. At the basis of these applications are the properties of 
the potential governing the gravitational and centrifugal forces operating in the system. In 
particular, the shape of equipotential surfaces and the existence of stationary points (where 
the net force exerted on a test particle vanishes) play a central role in our understanding of 
the evolution of binary star and planetary systems. 

More often than not, applications of the Roche model are built on the assumption that 
the orbit of the system is circular and that the system components are rotating synchronously 
with the orbital motion. Under this assumption, the components can be treated as static with 
respect to a co-rotating frame of reference and their shapes are determined by equipotential 
surfaces of the system. However, the assumption of circular orbits and synchronous rotation 
cannot always be justified, neither on observational nor on theoretical grounds. 

Observational support for non-circular orbits and non-synchronous rota tion is ve r y wel l 



established at present. Catalogs of eclipsing binaries as early as those by IShapleyl (119131 ) 



as the largest measured orbital eccentricity has since increased drastically (e.g. 


S 


ruve 


1950; 


Petrova & Orlov 


1999: Bildsten et al. 


1997: Heeediis et al. 2005: Pourbaix et al. 


2004 


) and 


now also includes binaries in the Small and Large Magellanic Clouds (e.g. 


Bavne et al. 


2002; 



Clausen et al.l 120031 ; iBayne et al.l 12004 ; iHilditch et al.l 120051 ) . Evidence for non-synchronous 
rotation of the components of close binary systems is equally abundant, and points to 



both sub- and supersynchronously rotating component stars fe.g. IStruvdll950l ; lLevatolll974 
Habets fc ZwaarJll989l : iMeibom etaDbood : Ivan Hamme fc Wilson! ll99oF l. 



Non-circular orbits are also observed in extrasolar planetary sys tems. Since the discov 



ery of the first extras olar planet orbiting a solar-type star, 51 Peg, (IMayor &: Quelozlll995 



Marcy fc Butler! Il995l ). the sample of known exoplanets has grown to more than 200, many 
of which orbit their host star with non-negligible orbital eccentricities. The most eccentric 
exoplanetary orbi t known to date is that of HD 20782b which has an eccentricity of 0.9 
( 1 Jones et al.ll2006l ). Constraints on exoplanetary rotation rates are yet to be inferred from 
observations. 



Theoretically, circularization and synchronization processes are driven by tidal inter- 
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actions between the components of binary or planetary systems^. Predicted circularization 
and synchronization time scales, however, depend strongly on the initial syst em parameters 



and theoretical uncertainties in the strength of tidal d issipation mechanisms (IMathieu et al. 
2004 ; iMeibom fc Mathieul 120051 ; iMeibom et al.ll2006l ). In addition, since synchronization 



tends to occur faster than circularization, components of binary and planetary systems tend 
to become synchronized with the orbital motion near periastron before the completion of 
circularization. The stellar and planetary rotation rates then deviate from synchronism at 
all other orbital phases. 

Such deviations from synchronism in gravitational two-body s ystems induce tim e-dependent 
oscillations in the atmospheres of the component stars or planets (IZahnlll97CH . Il975l ). There- 
fore these objects can no longer be treated as static with respect to a rotating frame of 
reference, unless the time scale of the oscillatory motions is sufficiently long compared to 
the dynamical ti me scale of th e star or planet. The validity of this approximation was 
first discussed by 
Kruszewski 1963; 



Limber! (119631) for the case of non-synchronous circular binaries (see also 



Savonijelll978l ). He concluded that as long as the rotational angular ve- 



locities of the stars do not deviate considerably from synchronicity, the components may be 
approximated as static and t heir shapes m ay be determined by the instantaneous equipo- 
tential surfaces of the binary. iLimberl (119631 ) also found the approximation to be valid in the 
limiting case of stars rotating with angular velocities close to the break-up angular velocity 
due to the dominance of the centrifugal force over the gravitational forces. 



Plaved (119581 ) was the first to abandon the assumption of synchronous rotation and 



study the effects of the asynchronicity on the shapes of the binary components and the posi- 
tion of the inner Lagrangian point L\ (the stationary point of the Roche potential located in 
between the two stars on the line connecting their centers of mass) in binaries with circular 
orbits. He found that the position of the L\ point tends to move closer to the mass center 
of a rapidly rotating binary component with increasing values of its rotational angular ve- 
locity. Correspondingly, the critical radius at which one component starts to transfer mass 
to its companion decr eases wi t h incr easing rotational angular velocity. Later investigations 
extended the work of iPlaved (119581 ) b y more carefully considering th e validity of the as- 



sumptions underlying the Roche model (lLimberlll963l ; Kruszewski! 1 19631) and consider i ng the 
effect of spin-orbit misalignment in circular binaries (lAvni fc Schiller! Il982l ). iLubowl (119791 ) 
also incorporated the effects of gas dynamics and heat transport in the surface layers of the 
binary components. 



For very close systems, orbital angul ar momentu m losses due to gravitational wave emission may also 
contribute to the circularization process |Peters!ll96j) . 
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Avnil (119761 ) further generalized the Roche model to account for eccentric binary orbits 



with the aim of modeling the light and radial-velocity variations observed in the e clipsing 



binary pulsar Vela X-l. His generalized binary potential was subsequently used by [Wilson 



( 119791 ) to discuss the computation of light and r adial- velocity curve s of eccentric binaries 



with non-synchronously rotating component stars. iRegos et al.l (120051 ). performed smoothed 
particle hydrodynamics calculations to validate the application of the Roche model to ec- 
centric binaries with component stars rotating close to synchronism with the orbital angular 
velocity at periastron. Part of their work, however, is flawed due to an incorrect term in the 
expression for the gravitational potential (see § [6] for details). 

Despite the vast observational support for the occurrence of eccentric orbits in binary 
and exoplanetary systems, a detailed account of the properties of the generalized Roche po- 
tential used to study these systems seems to be lacking from the literature. The geometry 
of the equipotential surfaces as well as the existence, location, and potential height of La- 
grangian points are necessary, for instance, in determining when mass transfer is initiated 
in non-synchronous, eccentric binaries along with the determining how mass and angular 
momentum flow during these mass-transfer phases. Consequently the theoretical study of 
component interactions in non-synchronous binaries of arbitrary eccentricity requires the 
understanding of this geometry and properties of the binary potential. Such interactions 
(e.g., mass exchange, mass and angular momentum loss) very often occur as binaries evolve 
and their modeling enters almost all population synthesis calculations of binaries of current 
interest (X-ray binaries, gamma-ray burst progenitors, binary compact objects, etc.). They 
are also relevant to the evolution of planetary systems along with the possibility of mass 
transfer and loss from planets in orbit around their host stars. The pursuit of theoretical 
understanding of these interaction processes in the context of broad binary evolution studies 
is what motivates the study presented here. 

The plan of the paper is as follows. In §m we introduce the Roche potential describing 
the motion of mass elements in non-synchronous, eccentric binary and planetary systems. 
In §|3l we discuss the time-dependence of the potential and investigate the conditions under 
which the potential varies sufficiently slowly to be considered quasi-static. Adopting these 
conditions, the existence and stability of stationary (Lagrangian) points of the potential are 
examined in §HJ In §51 we calculate the volume-equivalent radius o f the Roche Lo be as a 



function of the system parameters and provide a generalization of the lEggletonl (119831 ) fitting 
formula appropriate for non-synchronous, eccentric binary star and planetary systems. In 
§|f)l we examine the height of the potential along the line connecting the mass centers of the 
stellar or planetary binary components and discuss the implications for mass transfer and 
mass loss from the system. The final section, §3, is devoted to concluding remarks. 
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2. The Roche Potential for Non-Synchronous Eccentric Binaries 

We consider a binary system of starg^l with masses M\ and M2 orbiting one another in 
a fixed Keplerian orbit with period P rb, semi-major axis a, and eccentricity e. The distance 
D between the mass centers of the two stars and the Keplerian orbital angular velocity ujk 
of the stars are, at any time t, given by 

/ x «(l-e 2 ) 

D(t) = — -, 1 

V ; 1 + ecosz/' v ; 

and 

where 1/ is the true anomaly. 

We assume the first star (hereafter star 1) rotates uniformfyB about an axis perpendicular 
to the orbital plane with a fixed rotational angular velocity fli parallel to and directed in the 
same sense as the orbital angular velocity, ujk- We furthermore assume star 1 is sufficiently 
centrally condensed to be approximated by a point mass surrounded by a uniformly rotating 
zero-density envelope. We assume this envelope is convectively stable so that no bulk motions 
of the mass elements of star 1 occur other than those due to the star's rotation and the tidal 
force exerted by the companion. The companion star (hereafter star 2) is treated as a point 
mass. 



2 For brevity, we will simply refer to the system components as stars, bearing in mind that the analysis is 
valid for planetary systems as well. 



T he assumption of uniform rotat ion is made b ecause it is the simplest possible, though it is not without 



merit. 



Spruit fc Phinnev ( 1998h and Spruit (1999) argue that magnetic fields could establish uniform rota- 



tion. On the other hand, th ere are compelling case s of single and binary stars for wh ich differential rotation 
has been claimed (See, e.g., Popper fc Plavec 1976 : Domiciano de Souza et al. 2003^ . 
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Fig. 1. — Schematic representation of the position vectors used in the derivation of the 
equation of motion of a mass element (open circle) in a component of a non-synchronous 
eccentric binary. 
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The mass elements of star 1 are, at any given time, subjected to a sum of gravitational 
and rotational forces caused by the stars' mutual gravitational attraction, rotation, and 
orbital motion. We describe the motion of the mass elements of star 1 with respect to a 
Cartesian coordinate frame OXYZ with origin O located at the center of mass of star 1 
and Z-axis parallel to Qi. We furthermore let the frame rotate about the Z-axis with the 
rotational angular velocity of star 1. The equation of motion of a mass element of star 1 
with respect to OXYZ is then given by 

h = —VP - Wi - 2 fii x ?i, (3) 
P 

where fi is the position vector of the mass element with respect to the mass center of star 1, 
p and P are the mass density and pressure of the mass element, and 

y 1 = _ G ^_ G ^_W(X* + y*) + G^f (4) 
rd r 2 2 v ' JJ Z 



is the potential givin g rise to the gravitat ional and centrifugal forces acting on the mass ele- 



ment (for details, see lSepinsky et al.ll2007l . and references therein). Here, G is the Newtonian 
constant of gravitation, r*2 is the position vector of the considered mass element with respect 
to the mass center of star 2, and X and Y are the Cartesian coordinates of the mass element 
in the plane of the orbit. 

We note that Eq. (prj is valid only when the X-axis of the reference frame coincides with 
the line connecting the mass centers of the binary components. Since the initial orientation 
of the X- and F-axes in the orbital plane is arbitrary, the X-axis may be chosen to coincide 
with the line connecting the mass centers of the component stars at any given time without 
loss of generality. The definitions of the position vectors introduced in Eqs. ([3]) and (JH) are 
shown schematically in Fig. [TJ 

The potential V\ can be rewritten by expanding the gravitational attraction exerted by 
star 2 on the considered mass element of star 1 in terms of Legendre polynomials P n {x) as 

where R = fx — f 2 is the position vector of the mass center of star 2 with respect to the 
mass center of star 1, and x is the angle between R and f±. The first term in the expansion 
of — GMil\r<2\ is independent of the position of the considered mass element and thus does 
not contribute to the force exerted on the element. The second term in the expansion of 
— GMil\r^\ is identical to the last term in Eq. (j3j). Equation (j3J) is therefore equivalent to 

h = -^VP - VK - 2 fix x fi, (6) 
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where 



V* 



-G 



Mi 

In I 



\n 1 \ 2 {x 2 + Y 2 ) + w(n,t). 



(7) 



In the latter equation, W(fi,t) is the tide-generating potential defined as 



w(n,t) 




Pn (COS X) 



(.e.g. 



Zahnlll970l : IPolfliet fc Smeverslll99oh . 



Denoting the radius of star 1 by 1Z\ and restricting the tide-generating potential to its 
dominant n = 2 terms, it follows from dimensional analysis that 

2 



and 



where fL 



(l/2)\n 1 \ 2 (X 2 + Y 2 

w 




TZi 
D 



M. 
Mi 



(9) 



(10) 



GM 1 /\f 1 \ 

(GMx/TZfj 112 is the break-up an gular velocity of star 1. The centrifugal distor- 



tion thus becomes increasingly important when \Vti 
important when 1Z\ — > D and M 2 — > Mi. 



fL, while the tidal distortion becomes 



3. The quasi-static approximation 

For a binary with a circular orbit in which star 1 is rotating synchronously with the 
orbital motion, the potential V\ is independent of time and the mass elements of star 1 are 
at rest with respect to the co-rotating frame of reference. Equation ([3]) therefore reduces to 
the condition for hydrostatic equilibrium 

-VP = -W 1; (11) 
P 

so that surfaces of constant pressure and density coincide with surfaces of constant V\. The 
possible shapes of star 1 are thus determined by the equipotential surfaces of the binary. 

For a binary with an eccentric orbit or non-synchronously rotating component stars, 
Vi is a periodic function of time. At any given time, the mass elements of star 1 therefore 
undergo tidally induced oscillations on a time scale r tide determined by the difference between 
the instantaneous orbital angular velocity uJK{t) and the rotational angular velocity \Vti\: 

2tt . . 

Ttide = 1 m id if (12) 

\LO K (t) - &2i|| 
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From dimensional analysis, it follows that 



l r i| 



'dyn 



|(VP)/ P | |VVi| Tj 



(13) 



tide 



|Oi X r*i| X fi\ \Qi\ Tdyn 



(14) 



|(VP)/p| |Wi| r tide 

where Td yn = (GMi/lZf ) -1 / 2 is the dynamical time scale of star 1. Hence, when Td yn <C r t id e , 
the r*i and fii x r\ terms in Eq. ([3]) are negligible compared to the other terms, so that 
Eq. ([3]) approximately reduces to Eq. (TTTT) . In this limiting case, the motions of the mass 
elements are so slow that, at each instant, star 1 can be considered to be quasi-static with 
respect to the rotating frame of reference OXYZ. The instantaneo us shape of the star can 
then be approximated by the instantaneou s surfaces of constant V\ . iLimberl (119631 ) referred 
to this as the first approximation (see also Savonije 1978 ; Wilson 19791 ). 



The condition that Td yn <C r t id c can be expressed in terms of the orbital elements and 
properties of star 1 by writing r ti d e as 



Ttido 



where 



2tt 

Up 



Up 



1 + e cos v 



2tt (l + e) 1/2 
Porb (1 - e) 3/2 



/ 



(15) 



(16) 



is the orbital angular velocity at periastron, and / = |f2i|/Vp is the rotational angular 
velocity of star 1 in units of up. It follows that Td yn r t id c when 



P 



orb 



Tdyn 



> a(e, /, u), 



(17) 



where a is a function of the orbital eccentricity e, the ratio of rotational to orbital angular 
velocity at periastron /, and the true anomaly v: 



a(e, f, v) 



'l-e 



^2 
3/2 



1 + e cos v 
1 + e 



f 



18) 
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Fig. 2. — Variations of a(e, f, v) as a function of e, for v = and 0.8 < / < 1.2. Due to the 
symmetry of a(e, f, 0) about / = 1, the curves associated with / = 0.80, 0.90, 0.95, 0.99 are 
identical to those associated with / = 1.20, 1.10, 1.05, 1.01, respectively. For a(e, f, 0) < 1 
(below the dotted line), the quasi-static approximation is valid when P orb ^> r dyn . The inset 
shows the variation of 7d yn as a function of st ellar mass for stars on the ZAMS according to 
the mass-radius relation of iTout et al.l (119961 ). 



When the companion is at the periastron of the relative orbit, a(e, f, v) reduces to 



Hence, when / = 1 and the binary component stars are located at the periastron of their 
relative orbit, star 1 can be safely assumed to be quasi-static for all values of the orbital 
eccentricity e < 1. However, once / deviates from unity, the rapid increase of ujp with 
increasing orbital eccentricity causes a(e, f, 0) to also increase rapidly with increasing orbital 
eccentricity. This is illustrated in more detail in the main panel of Fig. [21 For the considered 
range of /-values, the function a(e, f, 0) is smaller than 1 for e < 0.6. Hence, in this 
parameter range, the quasi-static approximation is valid at the periastron of the binary 
orbit if P orb r dyn . As is illustrated by the inset of Fig. [2J this condition is easily satisfied 
for zero-age main sequence (ZAMS) stars of mass Mi < 10 M , provided the orbital period 
is longer than approximately 10 hr. For Mi < 1 M this condition may be relaxed even 
further. Due to the nonlinear dependence of a(e, f, v) on z/, the constraint on P or b may be 
more or less restrictive away from the periastron of the relative orbit. 



In the case of a binary with a circular orbit and synchronously rotating component 
stars, there are five points at which the gravitational forces exactly balance the centrifugal 
forces. A mass element with zero velocity experiences no accelerations at these points and 
thus remains stationary with respect to the co-rotating frame of reference OXYZ. Three of 
these so-called Lagrangian points lie on the line connecting the mass centers of the binary 
components, and two lie on the tips of equilateral triangles whose bases coincide with the 
line connecting the mass centers of the component stars. 

For binaries with eccentric orbits and/or non-synchronously rotating component stars, 
the tide-generating potential introduces a time-dependence in the potential V\ and thus in the 
position of the stationary points (assuming these points still exist). To examine the existence 
and properties of these points, we adopt the quasi-static approximation, Td yn <C r t id e , and 
introduce the following dimensionless quantities: 




(19) 



4. Lagrangian Points 



4.1. 



Existence and Location 
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From Eq. (j3j), it follows that 



V D = X D -^-- 



where 



r D (r 2 D - 2X D + 1) 1/2 
1 -{X 2 D +YZ)(l + q)A(f,e,u), (20) 



A(f,e,u)= f} 1 + e) * (21) 
(1 + e cos v) 



This factor groups all dependencies of the potential on the degree of asynchronism, the 
orbital eccentricity, and the mean anomaly. In the particular case of a binary with a circular 
orbit and synchronously rotating component stars, A = 1. 

The dependence of A(f, e, v) on e and v is illustrated in Fig. El for / = 1. Since A(f, e, v) 
is linearly proportional to the square of the rotational angular velocity in units of the orbital 
angular velocity at periastron, the curves are easily rescaled for different values of /. For 
a given value of u, A(f, e, v) increases with increasing e, while for a given e, A(f, e, v) 
reaches a maximum at apastron and decreases towards periastron. A(f, e, v) is furthermore 
most sensitive to the value of the orbital eccentricity near v = tt, where it increases from 
A(f, e, v) — 1 for e = to A(f, e, v) w 10 4 for e = 0.9. Moreover, at periastron, A(f, e, v) < 2 
for all e < 1 and / < 1. 
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Fig. 3. — Variations of A(f, e, v) as a function e and v, for / = 1. Variations of A(f, e, v) 
for other values of / are obtained by rescaling the plotted curves by a factor of f 2 . 
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The stationary points of the binary are solutions of the vector equation VVd = with 
respective X, Y, and Z components 



X 



D 



q 



D 



+ 



(r 2 D - 2X D + 1) 
1 

(r 2 D - 2X D + 1) 



3/2 



3/2 



- (l + q)A(f,e,v) 



+ 1 = 0, 



(22) 



Yi 



D 



± + 

3 



r 



D 



(r 2 D - 2X D + 1) 



3/2 



- (l + q)A(f,e,u) 



J D 



D 



2X 



D 



1 



3/2 



0. 



(23) 



(24) 



where we have used the fact that r 2 D — 2Xd + 1 = (Jo — l) 2 + Y^ + is a positive quantity 
Since both terms inside the square brackets in Eq. (1241) are positive, the equation can only be 
satisfied by setting Zd = 0. Thus, any and all stationary points necessarily lie in the plane 
of the orbit, as in the case for circular binaries with synchronously rotating component stars. 
We note in passing that Lagrangian points outside the orbital plane have been shown to 
exist when the spin axis of star 1 is inclined with respect to the orbital angu lar momentum 
axis of the binary (see, e.g., lAvni fc Schiller! Il982l ; iMatese fc Whitmird 119831 ). 



Letting Z D = 0, we solve Eqs. (|22|) and 
two possible cases: Y D = and Y D ^ 0. 



for Xd and Yd by distinguishing between 



(%) Y D = 0. Setting Y D = Z D = reduces Eq. ([22]) to 



X 



q 



D 



x D \ 3 



+ 



X 



D 



\x 



-X D (l + q)A(f,e,u) + l = 0. 



(25) 



D 



Leaving aside the trivial solution Xd = 0, this equation has three real solutions which can 
only be determined numerically and which correspond to three co-linear stationary points 
on the line connecting the mass centers of the binary components. In accordance with 
the nomenclature for circular synchronized binaries, we will refer to these points as the 
Lagrangian points L 1; L 2 , and L 3 . The L\ point is the stationary point located between the 
two component stars (0 < Xd < 1), the L 2 point is the stationary point located opposite 
star 2 as seen from star 1 (Xd > 1), and the L 3 point is the stationary point located opposite 
star 1 as seen from star 2 (Xd < 0). 
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Fig. 4. — The position the co- linear Lagrange points L 1? L 2 , and L 3 along the line connecting 
the mass centers of the binary components as a function of the logarithm of the mass ratio, 
q, for different values of A. 
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The location of each of the three co-linear Lagrangian points on the X-axis of the 
OXYZ frame is shown in Fig. H] as a function of q and A. From the middle panel of Fig. H] 
we can see that, for a given value of A, Xn increases with increasing values of q. For a given 
value of q, on the other hand, X L1 decreases with increasing values of A. Additionally, X L1 
asymptotes to a constant value for both large and small q. This behavior can be understood 
from Eq. (|25|) . which, for small q, reduces to 

X L1 [A (X L1 - l) 2 - (X L1 - 1) + 1] « 0. (26) 

The only real solution to this equation for which < Xn < 1 is Xli ~ 0. Thus, in the limit 
of small q, X L1 — > 0. For large q, Eq. (1231) yields 

X L1 « A- 1/3 . (27) 

This approximation breaks down for A < 1, in which case Xn — > 1 for large values of q. This 
change in the behavior of Xn for large values of q is due to the fact that the gravitational 
force exerted by star 2 is no longer negligible as Xn — > 1. 

For L 2 , we can see from the top panel of Fig. H]that, for a given value of A, X L2 decreases 
with increasing values of q, and, for a given value of q, Xi2 decreases with increasing values 
of A. For small and large values of q, Xi2 furthermore asymptotes to a constant value 
depending on the magnitude of A. This may be understood from Eq. (125]) which, for small 
q, reduces to 

XI ~ 2 -^Xl 2 + *±1 Xl2 -^0. (28) 

This equation admits of one real solution for X^ 2 which is a function of A. The analytical 
expression for the solution is, however, rather intricate and does not contribute further to 
understanding the behavior of X^ 2 for small values of q other than confirming the dependence 
of the asymptotic value of Xfi on A. For large values of q, Eq. (|25|) yields 

X L2 « A- 1/3 . (29) 

This approximation breaks down for .4. > 1 for the same reason that Eq. (12 7p breaks down 
for A < 1, i.e., the gravitational force exerted by star 2 is no longer negligible when X^ 2 — > 1. 
Instead Xl 2 « 1 for large values of q and A > 1. 

In the case of L 3 , we can see from the bottom panel of Fig. [4] that, for a given value of 
A, X L3 decreases with increasing values of q, while, for a given value of q, X L3 increases with 
increasing values of A. In addition, asymptotes to zero for small values of q, and to a 
constant depending on A for large values of q. This may again be understood from Eq. (|25l) 
which, for small values of q, reduces to 

X L3 [AX 2 L3 -(2A+l)X L3 + A + 2] wO. (30) 
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The only real solution to this equation with Xl3 < is Xls ~ 0. For large values of q, on 
the other hand, Eq. (j25|) yields 



X 



L3 



-A~ 1/3 . 



(31) 



for all values of A. 



(ii) Yd 7^ 0. For Yd ^ and Zd = 0, Eqs. fl22|) and ff23l yield two non-linear, algebraic 
equations: 

' 1 = 0, (32) 



(rl - 2X D + if 2 



1 



(l + q)A(f,e,u)=0. 



r 3 D (rl - 2X D + if 2 
These equations have two real, analytical solutions for Xp and Yd given by 

-i 2/3 



(33) 



X 



D 



l + q)A(f,e, v)-l 
Y D = ±y/X D (2-X D ). 



(34) 
(35) 



The solutions correspond to two stationary points located on the tips of triangles whose 
bases coincide with the line connecting the mass centers of the binary component stars. 
In accordance with the nomenclature for circular synchronized binaries, we refer to these 
points as the triangular Lagrangian points L4 and L5. The L4 point is the stationary point 
associated with Yd > 0, while the L 5 point is the stationary point associated with Yd < 0. 
Real solutions to Eqs. and (135]) furthermore only exist when 



q > 



A 



A-l/l 



(36) 



This condition is always satisfied when A > 1. In the limiting case where e = and 
/ = 1, Eqs. ( 1341 and ( l35i) reduce to the usual formulae for the position of the L 4 and L 5 
points in the case of a ci rcular, synchronized binary, i.e. Xd = 1/2 and Yd = ±v / 3/2 (e.g. 
Murray k Dermottl 120001 1 . 
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Fig. 5. — The Xo and Yd coordinates of the leading triangular Lagrangian point, L 4 , as a 
function of q for different values of A. The dotted line terminates near q = 0.69 below which 
L 4 does not exist for A = 0.64. The X D and Y D coordinates of the trailing Lagrangian point, 
L5, is identical to that of L 4 , except that Yd has the opposite sign for L 5 as for L 4 . 
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The dependence of the Xp and Yd coordinates of L4 on q and A is illustrated in the 
bottom and top panels of Fig. [51 respectively. Since X L5 = X Li and Y L5 = —Y L4: , the 
dependence of the position of L 5 on q and A is similar to that of L 4 . For a given value of 
A > 1, both and increase with increasing values of q. The coordinates furthermore 
asymptote to constant values for both small and large values of q. The asymptotic behavior 
can be understood from Eq. which reduces to 



for large values of q. Thus, in the limit of small q, Xn — > when A > 1. For A < 1, the L4 
point no longer exists for all values of q. E.g., when A = 0.64 (the dotted line in Fig. [5]), the 
L 4 point only exists for q > 0.69 (see Eq. [36]). However, when the L 4 point exists, X L4 still 
asymptotes to the value given by Eq. (1381) for large values of q. Finally, for a given value of 
q, X Li and Y Li decrease with increasing values of A. In addition, it is interesting to note 
that for A > 1, L 4 and L5 are always located inside the binary orbit, while for A < 1, L4 
and L5 are located outside the binary orbit. 




(37) 



for small values of q, and 



0.5 A~ 2 I Z 



(38) 



4.2. Stability 
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Fig. 6. — Stability regions (gray-shaded) for the Lagrangian points L 2 (left) and L 4 , L 5 
(right) in the (q, ^4.)-plane. The L 2 stability region remains always at log A < 0, in agreement 
with the known results for circular, synchronous binaries (log A = 0) for which L 2 is always 
unstable. The stability regions for L4 and L5 are identical to each other since the stability 
conditions given by Eqs. fl39l) - (HTl) are dependent only on the square of the Y D coordinate. 
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Following the methodology described by iMurray fc Dermott J2000h (Chapter 3.7) for 



circular, synchronous binaries, we examine the stability of the Lagrangian points in eccen- 
tric, non-synchronous systems. We employ a linear stability analysis where we perturb and 
linearize the equation governing the motion of a mass element in the binary potential and 
follow the behavior of an element given a small displacement from one of the Lagrangian 
points. If the solution to the perturbed equation of motion is oscillatory, the Lagrange 
point is stable to small perturbations. All other solutions imply instability. For stability, 
then, it follows that all three of the following conditions on the second derivatives of the 
dimensionless potential Vd must be satisfied. These are 

B 2 > AC, (39) 



VB 2 - AC < B, (40) 
VB 2 - AC > -B, (41) 



where 



B = AA(1 + q) + V xx + V YY , (42) 
C = V XX V YY -V XY , (43) 

Vxx = d 2 VD/dXl), Vyy = d 2 Vu/dY^, and Vxy = d 2 Vo/dX D dY D are the second derivatives 
of the potential Vd with respect to Xd and Yd, and the square roots in Eqs. (1401) and (14TI) 
are assumed to be positive. Examining these conditions, we find that L\ and L 3 are always 
unstable for any combination of binary parameters, while L 2 , L4, and L 5 can be stable as 
well as unstable depending on the values of q and A. The regions in the (q, ,A)-plane where 
L2, £4, and L 5 are stable and are shown by the gray-shaded areas in Fig. [6j 

For q < 0.1, the L2 point is found to be stable for a narrow range of A- values very near, 
but always less than A — 1, while for q > 10, a broad stable region is found at A < 0.1 (Fig.0 
left panel). In the latter regime, the potential is very similar to that of a non-rotating single 
star, except for a small contribution due to the gravitational potential of star 2. Even though 
this contribution is small, it is significant enough to be the sole cause of this stable region. 
This is further supported by the fact that, in this regime, the potential at L 3 is identical 
to that at L 2 except for this small contribution, and yet L 3 is never stable throughout the 
considered parameter space. Thus, the second body, however small, is the sole cause for 
the stable regions of L 2 . We n ote that we do not find an y stable points for log A = 0, in 



agreement with the analysis by IMurray fc Dermottl (120001 ) for circular, synchronous binary 
systems. 

The stability region for the triangular Lagrangian points, L4 and L5, covers a wide range 
of q and A values which, for q < 0.1, converges into a narrow stability region near A — 1 
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(Fig.O right panel). Solutions at exactly ,4=1 exi s t only for logg < — 1.4 and logg > 1.4, in 
agreement with the analysis of iMurray fc Dermottl (120001 ) for circular, synchronized binaries 
(see their Eq. 3.145). It is apparent, that the L 4 and L 5 stability regions are much larger for 
eccentric, non-synchronous systems than for circular, synchronous systems. The existence of 
these stable points is important to study the trapping of test particles su ch as gas attempting 



to escape the system or even Trojan asteroids in planetary systems (see lFord fc Gaudill2006 



for a discussion on Trojan asteroids around extrasolar planets). It should also be noted 
that Lagrangian points within the stable regions shown in Fig. O need not remain stable 
throughout the binary orbit due to the dependence of A on the mean anomaly (see Eq 
and Fig. ED. 



5. The Volume-Equivalent Roche-Lobe Radius 

Under the quasi-static approximation, the maximum size of a star in a close binary is 
determined by the equipotential surface connected to the inner Lagrangian point L\. This 
surface is commonly referred to as the Roche lobe. Once a star fills its Roche lobe, any 
further expansion of the star or shrinkage of the lobe results in mass loss from the star. For 
a circular binary with synchronously rotating component stars, the shape and volume of the 
Roche lobe depend solely on the mass ratio, q, of the system. For eccentric binaries with 
non-synchronous component stars, the shape and volume of the Roche lobe also depend 
upon the eccentricity, the true anomaly, and the degree of asynchronism. As before, we 
group these dependencies into the parameter A(f, e, v) (see Eq. |2"TT) . 
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Fig. 7. — Instantaneous equipotential surfaces of Vd given by Eq. (!20l) in the Z = plane for 
different mass ratios q = Mi/M 2 and values of the parameter A(f, e, v). The contours shown 
correspond to the equipotential surfaces passing through the co-linear Lagrangian points L%, 
L/2, and L3. The locations of the triangular Lagrangian points L4 and L5, if they exist, are 
marked with crosses. Star 1 is located at (X D , Y D ) = (0, 0), and star 2 at (X D , Y D ) = (1, 0). 
The A = 1.0 panels correspond to circular binaries with synchronously rotating component 
stars. For A = 0.64, L 4 and L 5 only exist for q > 0.69 (see Eq. l36l) . 
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The variations in the shape and size of the binary equipotential surfaces are shown in 
Fig. [7J as functions of q and A. The contours shown correspond to the cross-sections of 
the equipotential surfaces passing through the co-linear Lagrangian points with the Zd = 
plane. The locations of the triangular Lagrangian points, if they exist, are marked with 
crosses. For large mass ratios, the equipotential surfaces resemble that of a rotating single 
star regardless of the value of A. For smaller mass ratios, the structure of the equipotential 
surfaces becomes more complex and the detailed shapes become more sensitive to the value of 
A. A particularly interesting observation is that for A > 1 the equipotential surface passing 
through L\ may "open up" and no longer enclose star 2@. This geometry may facilitate mass 
loss from the system through L 2 when star 1 fills its Roche lobe and transfers mass to star 2 
through L\. We will discuss this in more detail in the next section. For values of A > 10, 
the centrifugal term dominates the potential Vd, so that the equipotential surfaces resemble 
those of a rotating single star regardless of the mass ratio q. 



The equipotential surfaces shown in Fig. [7J differ considerably from those of iRegos et al. 



(120051 ) . In particular, the orbital parameters adopted in Figs. 3 and 4 of that paper yield 
^.-values of 1.4 and 1.8, respectively, which may be compared to the A = 1.5 and A = 2.0 
panels for q = . 5 in F ig. [7J The di fferences may be a ttributed to the erroneous use of the 



Pringle &: Wadd (119851 ) potential in IRegos et al.l (120051 ) which is only valid for binaries with 
circular orbits. 

In order to determine whether or not a star fills its critical Roche lobe, a full three- 
dimensional treatment of stellar and binary evolution is, in principle, required. Such a treat- 
ment is, however, prohibitively computationally expensive. Roche-lobe overflow is therefore 
usually studied using one-dimensional stellar and binary evolution codes in which the Roche 
lobe of a star is approximated by a sphere of equivalent volume. The star of radius IZi is 
then assumed to overflow its Roche lobe and transfer mass to its companion when TZi > TZi, 
where TZl is the radius of the sphere with volume equal to that of the Roche lobe. 



4 Note that this "opening up" of the Roc he potential has a lso b een seen in studies of radiation effects on 
the shape of the equipotential surfaces. See Schuerman ( 19721 ) and Vanbeveren (1977). 
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Fig. 8. — The volume-equivalent Roche lobe radius, TZl in units of the instantaneous distance 
D(t) between the binary components as a function of log q, for a range of A- values. Small 
fluctuations in the curves are due to the finite errors in the Monte Carlo integration method 
used to calculate the volume of the Roche lobe. 



-26- 



In Fig. [HI we show the volume-equivalent Roche lobe radius in units of the instantaneous 
distance between the binary components calculated via Monte Carlo volume integration as 
a function of the mass ratio, q, for a ra nge of A- values. W e note that our calculation of 



the Roche lobe radius agrees with that of IVanbeverenl (119771 ) to within a few percent in the 



regime in which our parameter spaces intersect (circular, non-synchronous orbits where the 
effects of radiation pressure have been neglected). Comparison of Fig. [8] with Fig. H] shows 
that TZi as a function of q is similar in shape to L\ as a function of q, which is expected 
since it is the potential at L\ which defines the size and shape of the Roche lobe. For small 
values of q, TZl asymptotes to zero independent of A because the location of L\ asymptotes 
to the position of star 1. For large values of q, TZl asymptotes to an ^4-dependent value. 

A simple and accurate fitting formula for the volume-equivalent Roche Lobe radius of 
a s tar in a c i rcular binary with synchronously rotating component stars has been provided 



bv lEggletonl (119831 ): 



T7 49 a 2/3 

L ' circ 0.6g 2 / 3 + ln(l + gV3)' ^ 

where a is the radius of the circular orbit. For lack of a better treatment, this formula is 
often extrapolated to binaries with eccentric orbits. In particular, Roche lobe radii of stars 
at the periastron of a binary orbit are commonly approximated by 

^^(l-e) 06g2 ; + 4 g 3 +gl/3) , (45) 

where a(l — e) is the periastron distance of the binary orbit. We can further extend this 
generalization and apply it to arbitrary orbital phases by replacing the periastron distance 
between the stars with the current orbital separation, D(t). Thus, 

7? Egg D(t) ( A Q\ 

Kl ^ W 0.6g 2 /3 + i n (i + g i/3)- ^) 
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Fig. 9. — The ratio of the volume-equivalent Roche Lobe radius at periastron obtained by 
numerically calculating the volume of the Roche lobe to the volume-equivalent Roche Lobe 
radius obtained from Eq. (|46p . Small fluctuations in the curves are due to the finite errors 
in the Monte Carlo integration method used to calculate the volume of the Roche lobe. 
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In Fig. [9J we show the ratio of the volume-equivalent Roche Lobe radius at periastron 
calculated through a Monte Carlo integration method to the volume-equivalent Roche Lobe 
radius determined by means of Eq. ()46[) . For A ~ 1, we find Eq. (|46j) to be accurate to within 
a few percent over a wide range of g-values, confirming previous statemen ts on the accuracy 



of the equation b ased on smoothed particle hydrodynamics calculations (IFaber et al.l 12005 



Regos et al.l 12005 ). For a given value of A, the ratio TZl/TZ^ is closest to unity near q = 1 
and indicates better agreement between IZl and 7Z^ ee for q < 1 than for q > 1, especially 
for large values of A (by almost a factor of 2 as shown in Fig. [9]). For a given value of q, the 
agreement between TZl and 7Z^ sg furthermore becomes progressively worse as A differs from 
unity. 

Hence, a star rotating synchronously with the orbital angular velocity at periastron in 
an eccentric binary (so that 1 < A < 2 for e < 0.9) has a smaller volume-equivalent Roche 
Lobe radius at periastron than a star rotating synchronously in a circular orbit with radius 
equal to the periastron distance of the eccentric binary. As shown in Fig. [9], this decrease in 
IZl may be substantially smaller than the (1 — e) factor introduced in Eq. (145]) . Consequently, 
the maximum size of star 1 is smaller in an eccentric binary and mass transfer from star 1 
to star 2 may begin earlier than in a circular binary with the same instantaneous distance 
between the component stars. 

For ease of use, we provide a fitting formula for TZl/TZ^ as a function of q and A. In 
order to keep the accuracy of the fit better than 1%, it is necessary to divide the parameter 
space into six different regimes. The formulae below are accurate to better than 1% for 
—8 < logg < 8 and —2 < log A < 4, and, in most cases, the accuracy is better than 0.5%. 



[logg > 0;\ogA < -0.1] = 

1.226 - 0.2L4 - 0.15(1 - A) exp[(0.254 - 0.3)(logg) L55 ] (47) 



TZ E L gs 



[\ogq < 0;log^l < -0.1] = 

1 + 0.11(1 -A)- 0.05(1 - A) exp[-(0.5(l + A) + logg) 2 ] (48) 



[log? < 0; -0.1 > log .4 < 0.2] = 9o (A) + gi(A) logg + g 2 {A){logqf 



(49) 
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g (A) = 0.9978 -0. 1229 log „4-0.1273(log„4) 2 
gi (A) = 0.001 + 0.02556 log A 
g 2 (A) = 0.0004 + 0.0021 log A 



^ [log g > 0; -0.1 > log^l < 0.2] = h (A) + h(A) log? + h(A)(logq) 2 (50) 

h (A) = 1.0071 - 0.0907 log „4 - 0.0495(log„4) 2 
h^A) = -0.004 - 0.163 log „4-0.214(log„4) 2 
h 2 (A) = 0.00022 - 0.0108 log A - 0.02718(log„4) 2 



Egg 



[logg < 0; log„4 > 0.2] = i (A) + h(A) e*P [-i 2 (A)(logq + i 3 (A)) 2 ] (51) 

6.3014(log^) L3643 



i (A) = 
h(A) = 

12(A) = 



exp[2.3644(log„4) - 70748 ] - 1.4413 exp[-0.0000184(log^)- 4 - 5693 ] 
log A 

0.0015 exp[8.84(log A) - 282 ] + 15.78 
1 + 0.036 exp[8.01(log^) a879 ] 



0.105 exp[7.91(log^) - 879 ] 

0.991 

13(A) = 



1.38 exp[-0.035(log^)°- 76 ] +23.0 exp[-2.89(log„4) a 



701 



[log g > 0;log.4> 0.2] =j (A) + ji(A)exp [-j 2 (A) (log q)^} (52) 

. 1.895(log^) a837 

30 ~ exp[1.636(logA) - 789 ] - 1 
. 4.3(log^)°- 98 

31 ~ exp[2.5(log„4) - 66 ] +4.7 

. _ 1 

32 ~ 8.8exp[-2.95(log„4) a76 ] + 1.64exp[-0.03(log„4)°- 76 ] 
j 3 = 0.256 exp[-1.33(log^l) 2 - 9 ](5.5exp[1.33(log^l) 2 - 9 ] + 1) 
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6. Systemic mass loss 



As mentioned in the previous section, once star 1 fills its Roche lobe, any further stellar 
expansion or orbital contraction may cause the star to begin losing mass through L\. In the 
standard case of a circular synchronous binary, any matter flowing from star 1 through Li 
will, in the absence of other non-gravitational or rotational effects, be captured by star 2, 
either forming a disk of material or falling directly onto its surface. In the case of an eccentric 
and/or non-synchronous binary, the geometry of the equipotential surfaces is such that this 
is not necessarily the case. 
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Fig. 10. — Height of the dimensionless potential Vd at periastron along the ^Tp-axis (Yd = 
Zd — 0) for different values of the mass ratio, orbital eccentricity, and ratio of rotational to 
orbital angular velocity at periastron. Different line types in each panel correspond to the 
orbital eccentricities e = (solid lines), e = 0.3 (dotted lines), e = 0.6 (short-dashed lines), 
and e = 0.9 (long-dashed lines). In each panel, the maxima in the curves from left to right 
correspond to the L 3 , L 1; and L 2 point. Star 1 is located at X D = and star 2 at X D = 1. 
For ease of comparison, the curves have been normalized to the height of the potential at 

u. 
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In particular, we have seen from Fig. [7] that for significantly large A, the equipotential 
surface passing through L\ no longer encloses star 2. Matter transferred from star 1 to star 2 
through the inner Lagrangian point Li may therefore be lost from the system. Whether or 
not this happens depends on the particle trajectories in the mass-transfer stream which we 
will address in detail in a forthcoming paper. For the purpose of this paper, we restrict 
ourselves to comparing the height of the potential at the three co-linear Lagrangian points 
as an indicator of possible mass loss scenarios in eccentric binaries. 

Since the first occurrence of mass transfer in eccentric binaries is expected to take place 
when the stars are closest to each other, we focus our discussion on the height of the potential 
at the periastron of the binary orbit. In addition, the graphical representation of the height 
of the potential turns out to be clearer in terms of the mass ratio, orbital eccentricity, and 
ratio of rotational to orbital angular velocity at periastron than in terms of the parameter 
A. Hence, the height of the potential Vd at periastron along the line connecting the mass 
centers of the component stars (the X^-axis) is shown in Fig. [10] as a function of q, e, and /. 
To facilitate the comparison of the relative height of the potential at the co-linear Lagrangian 
points, the potential is furthermore normalized to its value at the inner Lagrangian point 
L x . 

For the considered ranges of mass ratios, eccentricities, and rotation rates, the relative 
height of the potential at L2 depends most sensitively on e and /, while the relative height 
of the potential at L 3 is largely constant over the considered range of parameters. For a 
constant / and q, an increasing eccentricity tends to decrease the relative height of the 
potential. The decrease is stronger for larger /-values and smaller g-values. 

For the considered parameter ranges, the shape of the normalized potential near L\ is 
furthermore highly insensitive to the orbital parameters, while the height of the potential, 
as well as the position of the peak, can change for both L 2 and L 3 . Most interestingly, the 
potential at L2 may be higher than, equal to, or lower than the potential at the L\ point. 
Thus, matter flowing from star 1 through the L\ point will not necessarily be captured by 
star 2, and may be energetic enough to escape the system. As such, it is possible that 
mass transfer through L\ can be significantly non-conservative due to the geometry of the 
equipotential surfaces in eccentric binaries, but a calculation of the amount of matter lost 
from the system will require detailed particle trajectories, which is beyond the scope of this 
paper. We also note that the value of the potential at L3 can, for certain values of the orbital 
parameters, be very near to the value of the potential at L\. In cases such as this, it may be 
possible for a small amount of matter also to be lost from star 1 through L 3 in addition to 
that lost through L\. However, contrary to iRegos et al.l (120051 ). we find the potential at L 3 



to be always higher than the potential at L\, so that mass loss through L 3 never takes place 



pri or to mass loss through Ljj. This di fference may be at tributed to the erroneous use of 
the IPrinde fc Wade! Jl985h potential bv lReeos etaD J2005h . 



7. Summary 

In this paper, we present a detailed investigation of the existence and properties of 
equipotential surfaces and Lagrangian points in eccentric binary star and planetary systems 
with non- synchronously rotating component stars or planets. The analysis is valid as long 
as the time scale of any tidally induced oscillations is long compared to the dynamical time 
scale of the component stars or planets. Once the foundation of the potential is laid, we (i) 
solve for the positions of the three co-linear and two triangular Lagrangian points, (ii) study 
the stability of the points, (iii) provide a semi-analytical fit for the general solution of the 
volume-equivalent Roche lobe radius that allows the determination of when mass transfer is 
initiated, and (iv) discuss the role of the potential geometry and structure in inducing the 
possibility of systemic mass loss. Throughout the analysis we compare our results to the 
well-studied case of synchronous binaries with circular orbits and quantify the differences, 
which turn out to be significant depending on the binary properties. 

Our study has been motivated by the need for characterizing the gravitational poten- 
tial for the general case of eccentric, non-synchronous binaries when considering interactions 
(mass transfer or mass and angular momentum loss) that occur in binaries before full syn- 
chronization and circularization has been achieved. The analysis presented here can be 
applied to both stellar and planetary systems for a wide range of problems. We plan to use 
the results to study the ballistic motion of test particles in the potential, assess the conditions 
under which mass transfer may be conservative or not, and eventually apply the results in 
the development of a theoretical framework that allows modeling of the orbital evolution of 
interacting, eccentric, non- synchronous binaries. 

We thank Chris Belczynski and Fred Rasio for useful discussions, as well as an anony- 
mous referee for the encouragement to complete our analytic fits for the Roche lobe radius.. 
This work is supported by a NASA GSRP Fellowship to JS, a Packard Foundation Fellowship 
in Science and Engineering and a NSF CAREER grant (AST-0449558) to VK. 



5 Note that 



Regos et al. 



(|2005l ) adopt a nomenclature for L2 and L3 that is opposite to ours. 
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